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Abstract. We examine the cosmic ray proton (CRp) population within the giant elliptical galaxy M 87 using the TeV y-ray 
detection of the HEGRA collaboration. In our scenario, the y-rays are produced by decaying pions which result from hadronic 
CRp interactions with thermal gas of the interstellar medium of M 87. By comparing the y-ray emission to upper limits from 
EGRET, we constrain the spectral index of the CRp population to ajfy < 2.275 within our scenario. Both the expected 
radial y-ray profile and the required amount of CRp support this hadronic scenario. The accompanying radio mini-halo of 
hadronically originating cosmic ray electrons is outshone by the synchrotron emission of the relativistic jet of M 87 by one 
order of magnitude. According to our predictions, the future GLAST mission should allow us to test this hadronic scenario. 
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1. Introduction 

The giant elliptical galaxy M 87 is an intensively studied ob- 
ject in our direct extragalactic vicinity situated at a distance 
of 17 Mpc dNeilsen & Tsvetano\l200ol) . The announcement of 
the TeV y-ray detection of M 87 at a 4-cr s ignificance level by 
the HEGRA collaboration ( Aharonian et al. 2003) using imag- 
ing atmospheric Cerenkov techniques was the first discovery of 
TeV y-rays from a radio galaxy with a jet whose axis forms a 
relatively large angle wi th the line of sight of roughly 30° - 35° 
feicknell & Begermanll99ol) . On the basis of the limited event 
statistics the detected emission is inconclusive whether it orig- 
inates from a point source or an extended source. Despite test- 
ing for burstlike behavior of M 87 no time variation of the TeV 
y-ray flux has been found. This detection provides the unique 
possibility for probing different y-ray emission scenarios and 
thus provides new astrophysical insight into high energy phe- 
nomena of this class of objects. 

In the literature, there are three different types of model 
predicting y-ray emission from objects like M 87: In the first 
scenario, the GeV/TeV y-ray emission is generated by the 
active galactic nucleus (AGN), and possibly relat ed to pro - 
cessed radiation of the relativistic outflow ( Der mer et al.l 19971) . 
Particularly, inverse Compton (IC) scattering of cosmic mi- 
crowave background photons off electrons within the jet which 
have been directly accelerated or reaccelerated as well as the 
Synchrotron Self Compton scenario could lead to y-ray emis- 
sion (Bai & Ledl200ll) . Secondly, dark matter annihilation or 
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decay processes could be another conceivable source of y- 
ray emission, such as the hypothetical neutralino annihila- 
tion JBaltz et al.ll2000h . Finally, hadronic cosmic ray proton 
(CRp) interactions with the thermal ambient gas would pro- 
duce pion decay induced y-rays as well as inverse Compton 
and sy nchrotron emiss ion by secondary cosmic ray electrons 
(CRe) jVestrandll983l . These processes are possible due to 
the long lifetimes of CRp comparable to the Hubble time 
JVolk et.alJll99r3l . long enough to diffuse away from the pro- 
duction site and to maintain their distribution throughout the 
cluster volume. Because of the strong dependence of this 
hadronic process on particle density, the giant radio galaxy 
M 87, located inside the central cooling flow region of the 
Virgo cluster, is expected t o be a major site of y-ray emission 
JPfrommer & Enfilinl2003l) . 



This work uses the hadronic scenario to model the result- 
ing y-ray emission. Thus it probes the CRp population by the 
recent TeV y-ray observations yielding either an upper limit or 
a detection on the CRp population, provided this scenario ap- 
plies. However, this approach only constrains the CRp within 
the central region of intracluster medium (ICM) of the Virgo 
cluster which is dominated by the interstellar medium (ISM) 
of the radio galaxy M 87. In the following, we use the term 
ICM for both. It should be emphasized that this hadronic sce- 
nario predicts stationary y-ray emission and will be ruled out if 
the emission is found to be time-variable (barring the existence 
of a second component). This, however, would result in even 
tighter constraints on the CRp population owing to the absence 
of inescapably accompanying y-ray emission. 
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2. y-ray emission from hadronic CRp interactions 

The differential number density distribution of a CRp popula- 
tion can be described by a power-law in momentum p p , 



/ p (r, p p ) dp p dV = HcRp(r) 



GeV J 



(Ppf_\ 
\ GeV/ 



dV, 



(1) 



where the normalization ncRp(f) is determined by different 
models of spatial distribution of the CRp population (see 
Sect.rj. 

If the CRp population within the cooling flow region had 
time to loo se energy by means of Coulomb interactions in 
the plasma ( Gould 19721) . the low energy part of the spectrum 
would be modified. This can be treated approximately by im- 
posing a lower momentum cutoff 

P*=fhMC*2X>(±)(^)Q*Vc-\ (2) 

where we inserted typical values for M 87. The kinetic energy 
density of such a CRp population is 



SCRpO) 



«CRpO") m p c 2 irripC 



2(a-l) \GeV 



(3) 



1 + p 2 - 1 



where SJa,b) denotes the incomplete beta-function, x = (1 + 
p 2 Y l , and p = £2^. 

* 1 ' r m p c 

The CRp interact hadronically with the thermal ambient gas 
and produce pions, provided their momentum exceeds the kine- 
matic threshold p^ = 0.78 GeV c _1 of the reaction. The neu- 
tral pions decay into y-rays while the charged pions decay into 
secondary electrons (and neutrinos). Only the CRp population 
above the kinematic threshold p t h T is visible through its decay 
products in y-rays and thus constrained by this work while its 
lower energy part can not in general be limited by considering 
hadronic interactions only. 

An analytic formula describing the omnidirectional (i.e. in- 
tegrated over 4 n solid angle) diffe rential y-ray source functio n 
resulting from 7r°-decay is given in Pfrommer & Enfilin ( 2003): 



,_ ncRp(f) 
q 7 (r, Ey) dE y dV - cr pp c n N (r) 2 a 

Uev 



(4) 
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dE y dV , 



where n^(r) is the target nucleon density in the ICM assuming 
primordial element composition. The formalism also includes 
the detailed physical processes at the threshold of pion produc- 
tion like the velocity distribution of CRp, momentum depen- 
dent inelastic CRp-p cross section, and kaon decay channels. 
The shape parameter 5 and the effective cross section cr pp de- 
pend on the spectral index of the y-ray spectrum according to 



6 = 0.14ff" L6 + 0.44 and 
o- pp = 32-(0.96 + e 4 - 4 - 24ff ) mbarn. 



(5) 
(6) 



Provided the CRp population has a power-law spectrum, 
the relation of the hadronic y-ray flux T y in different energy 



bands can easily be found using the analytic formulae for the 
integrated y-ray source density (Pfrommer & Enfilin 2003), 



T y {E\ < Ey < E 2 ) A 7 (E U E 2 ) 



T y (E 3 <E y < E A ) 
where A y {Ei,Ej) = 

and Xij = 
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Here we introduced the abbreviation = f(xi) - f(xi). 

This relation for hadronic y-ray fluxes is independent of any 
specific model of CRp spatial distribution as long as the same 
physical mechanism governs the CRp distribution in both en- 
ergy bands. 

Using the HEGRA y-ray flux for M 87 of TjjE > 
730 GeV) = 9.6 x 10 I3 y cm 2 s 1 dAharonian et alJl2003h . 
and requiring the expected y-ray flux above 100 MeV to be 
smaller than the EGR ET upper limit Td E > 100 MeV) = 
2.18x 10~ 8 y cirT 2 s~' (Reimer et al. 2003), we are able to con- 
strain the CRp spectral index to a < 2.275. For this calcula- 
tion, we assume a constant CRp spectral index a<? v extending 
from the GeV to TeV energy regime. In the case of steeper 
spectra in the TeV region, the CRp spectrum needs to be bent 
in a convex fashion or to exhibit a low energy cutoff in order to 
meet the requirement imposed by EGRET. 

3. Spatial distributions of CRp within the ICM 

In the following, we introduce three models for the spatial dis- 
tribution of CRp within the ICM. The origin of the CRp popu- 
lation is not specified in the first two models, but the CRp may 
be accelerated by shock waves of cluster mergers, accretion 
shocks, or result from supernova driven galactic winds. 

The isobaric model assumes that the average kinetic CRp 
energy density £cRp(f) is a constant fraction of the thermal en- 
ergy density SthW of the ICM 



£CRp(r) = XcRp e th (r) . 



(10) 



The thermal energy density s,h( r) is obtained from a sp herically 
symmetric temperature profile dPfrommer & Enfilinl2003l and 
a double g-mode l of electron dens ities adapted to X-ray obser- 
vations of M 87 (Matsushi ta et al.l20 02). 

The adiabatic model assumes the CRp population to be 
originally isobaric to the thermal population but to become adi- 
abatically compressed during the formation of the cooling flow 
without relaxing afterwards: 

fiCRpW = Xcrp e th (r) -> £cR P ( r ') = ^CR P ( r ') £ th(r') • (11) 

Adiabatic compression of CRp changes the scaling parameter 
XcRp to 



X CRp( r ) = X CRp 



n e (r) 



(<*+2)/3 



- ^CRp 



T \(»+2)/3 
' cluster 1 



, (12) 



where r c ] ustel denotes the electron temperature in the outer re- 
gion of Virgo. The last step assumes hydrostatic equilibrium of 
the gas during this transition. 
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Diffusion of CRp away from M 87: The relativistic 
plasma bubbles produced by M 87 likely contain relati vistic 
proto ns, which can partly escape into the thermal ICM ( Enfilin 
2003|. Most of the CRp that have been injected into the clus- 
ter center are either diffusively transported into the surrounding 
ICM or form relativistic bubbles which rise in t he gra vitational 
poten tial of the cluster due to buoyant forces JChurazov et alJ 
l200ll and references therein). 

Momentum dependent CRp diffusion in a turbulent mag- 
netic field with a Kolmogorov-type spectrum on small scales 
would result in spectral steepening and therefore would violate 
the limits on the spectral index aSX obtained in Sect. |2] pro- 
vided there is no sharp upper cutoff in the CRp spectrum. Thus 
we adopt for simplicity the scenario of passive advective trans- 
port of CRp in a turbulent flow with a diffusion coefficient k 
independent of momentum. The time-dependent CRp distribu- 
tion function reads for short (duration At) point-like injection 
with CRp injection rate Q(p p ) at time t = 



f P (r,p P ,t) 



g(p P )Af 

{AntKfl 2 6XP \ At K 



(13) 



In a quasi-stationary situation, which is a valid approximation 
for timescales longer than the typical CRp diffusion timescale 
in the case of a stationary or short-term intermittent CRp 
source, the integrated CRp distribution function is approxi- 
mately given by 



f p (r,p p )= f dt'f p (r,p p ,t') = 
Jo 



QiPy) 
4n tk 



erfc 



(14) 



where erfc(x) denotes the complementary error function which 
is responsible for the spatial cutoff at the characteristic diffu- 
sion scale /?diff = V6 1 k. While assuming a power-law shaped 
Q(p p ), the CRp distribution of eq. Q can be written within the 
framework set by this model as 



«CRp,o c / r \ 




(15) 



Following Pfrommer & Enfili^ J2003I) . the averaged CRp 
luminosity of M 87 can be estimated to be 



L C Rp 47rm p c 2 ncR P ,okpc (m p c 2xl 



2(a-l) 



GeV 



& 



a-2 3 



(16) 



Considering an energy dependent mean-free-path for diffusion 
would require fine-tuning of this model while simultaneously 
enlarging the accessible parameter space. While the resulting 
profiles should only be slightly affected by this change, this 
could possibly alter our conclusions concerning the normaliza- 
tion «cRp,o- However, this would only add further uncertainty 
to «cRp,o which already depends on two unknown parameters, 
the lifetime of the source t and the diffusion coefficient k. 

4. Modeled y-ray profiles 

The y-ray flux profiles f 7 (r ± ,E 7 > Zs t hr) are obtained by inte- 
grating the y-ray source function q y (r, E 7 ) of eqn. © above a 
threshold energy E$a, successively projecting and convolving 
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Fig. 1. a) Modeled y-ray surface flux profiles f 7 (r±) as function 
of impact parameter r ± in our three different models for the 
spatial distribution of the CRp population. They are normal- 
ized by comparing the integrated y-ray flux above 730 GeV to 
HEGRA data of M 87 within the innermost two data points, b) 
Comparison of detected to integrated y-ray flux T 7 {6 2 ) within 
the central aperture and the innermost annuli for different mod- 
els of spatial CRp distribution as well as two different widths of 
the PSF. The thick black lines correspond to cr — 0.05° whereas 
the thin grey lines are calculated for cr = 0.08°. The vertical 
dashed line separates the data from the noise level at a position 
corresponding to r ± = 37.5 kpc. 



the spherically symmetric profiles with the point spread func- 
tion (PSF) of HEGRA, PSF(rJ = exp[-r 2 J(2 o- 2 )]/(2n cr 2 ). 
Resolution studies based on observations of the Crab Nebula 
with HEGRA indicate a width of cr = 0.08°, assum ing a differ- 
ential spectral index of a — 2.7 (Dau m et alJll997l) . However, 
for flatter power-law spectra being preferred by our hadronic 
y-ray model (see Sect. the width will be smaller owing 
to increasing mean y-ray energy. This leads to an increase of 
y-ray induced particles of the air shower and therefore better 
quality of shower reconstruction according to a smaller relative 
Poissonian error. For a rough estimate, we rescaled the width 



4 



C. Pfrommer & T. A. Enfilin: Probing the Cosmic Ray Population of the Giant Elliptical Galaxy M 87 



of the PSF using the scaling of the mean y-ray energies above 
the instrumental threshold, yielding cr = 0.05° with a = 2.2. 

The line-of-sight integration was performed out to a radius 
of R m ax — 3 Mpc which corresponds to the characteristic dis- 
tance where the yS-model of electron densities is no longer ap- 
plicable due to accretion shocks of the cluster. The resulting 
y-ray profiles are shown in Fig. ^ The normalization of the 
surface fluxes depends on the assumed scaling between CRp 
and thermal energy density, which is fixed by comparing the 
integrated flux above 73 GeV to the innermos t two y-ray flux 
data points of HEGRA ( Aha ronian et alJl20 03) corresponding 
to 6 — 0.126° or r± = 37.5 kpc. Although there are distinct 
morphological differences visible in the three spatial CRp mod- 
els, the convolution with the PSF leads to very similar profiles 
for the expected HEGRA y-ray counts within the uncertain- 
ties. To demonstrate this, we compare the integrated y-ray flux 
Fy(E y > Etta-) f° r different annuli of equal solid angle elements 
centered on the source to the HEGRA data for the two different 
widths of the PSF discussed above (Fig.0. 

5. Consequences for the CRp population in M 87 

By employing the technique described in Sect. 0] we explore 
the consequences for the CRp scaling parameters Xqr p and 
"CRp.o in the particular models of CRp spatial distributions. The 
resulting values, shown in Tabled have been obtained using 
a PSF of width cr = 0.05°, however there are no significant 
changes in Xqr p for <x = 0.08°. The values of the CRp scal- 
ing parameter Xcr p inferred from M 87 are comparable to the 
one in our Galaxy, which is of order unity (Parker 1969). Since 
the HEGRA y-ray measurements probe only the central region 
of the Virgo cluster which is dominated by the elliptical ra- 
dio galaxy M 87, a composition of ISM and ICM is observed, 
potential ly mixed by c onvective motion within the cooling flow 
JChurazov et alJ200ll) . Therefore we expect Xqr p to be smaller 
than in our Galaxy, but significantly higher than upper limits 
obtai ned in nearby cooling flow clusters, which are less than 
20% (Pfrommer & EnRlinl200i . 

In the case of diffusion of CRp away from M 87, we are 
able to constrain the averaged CRp luminosity LcRp of the cen- 
tral AGN by assuming a plausible value for the diffusion co- 
efficient k. The inferred values are of the same order as in- 
stantan eous jet power estim ates of M 87, Lj et — 10 43 erg s _1 
jBicknell & Begelman 1996). Thus, we limit a combination of 
diffusion efficiency of CRp into the ambient thermal medium 
and average jet power by this approach. 

Because of the scaling behavior of Xcr v in the isobaric and 
adiabatic models, we quantify the influence of a lower cut- 
off Pmm on the population of CRp due to Coulomb interac- 
tions in the plasma by taking the ratio of CRp energy densities 
ecRp(Pmin) with and without lower cutoff (see eqn. 0). Such a 
cutoff yields lower values of Xcr v and therefore smaller contri- 
bution to the y-ray flux in the energy range of EGRET once the 
CRp momentum cutoff p m \ n exceeds the kinematic threshold 
Ahr = 0.78 GeVc" 1 of the hadronic interaction. Thus, cooling 
of the CRp population allows for steeper power-law distribu- 
tions. 



Table 1. Consequences for the CRp scaling parameter Xqr p 
and «cRp,o by comparing the integrated flux above 730 GeV to 
HEGRA data of the radio galaxy M 87 within the innermost 
two data points corresponding to 9 = 0.126°. The spatial dis- 
tribution of CRp is given by the isobaric, the adiabatic, and the 
diffusion model, respectively (see Sect. fJJ. In the first two cases 
the values are calculated for a CRp population with and without 
lower cutoff p m \ n while in the latter case two different lifetimes 
of the source have been considered. Note that the averaged CRp 
luminosity Lcr p scales with K29 = k/(10 29 cm 2 s _1 ). 



"GeV 


yisobaiic 
A CRp 


A CRp 


"cRp.o [cnr 3 ] 


^cr p [*29 erg/s] 




"mill 


GeVc-' : 


t = 1 Gyr k-_1 : 




2.1 

2.2 
2.27 


0.47 

0.65 
0.99 


0.31 
0.43 
0.64 


1.8 x 10~ 7 
5.0 x 10- 7 
1.0 x 10~ 6 


1.1 x 10 43 
1.6 x 10 43 
2.4 x 10 43 






2 GeVc- 1 : 


t - 3 Gyr : 




2.1 

2.2 
2.27 


0.42 
0.52 
0.73 


0.28 
0.34 
0.47 


1.2 x 10~ 7 
3.4 x 10~ 7 
6.8 x 10~ 7 


7.6 x 10 42 
1.1 x 10 43 
1.6 x 10 43 



6. Synchrotron emission by hadronic CRe 

Follow ing the formalism described in Pfro mmer & En filin 
( 2003J), we compute the synchrotron emission of CRe result- 
ing from hadronic CRp interactions. Integrating the expected 
radio surface brightness profiles over the solid angle element 
corresponding to the y-ray emission region and assuming mag- 
netic fields of the form B{r) = IOjuG [n e (r)/n e (0)] 05 , we ex- 
pect hadronic synchrotron fluxes F v = Fq [v/(1.4 GHz)] _q- ^ 2 , 
where Fy = 11 Jy and 16 Jy for a = 2.1 and 2.2. However, this 
hypothetical radio mini-halo is outshone by the synchrotron 
emission of the relativi stic jet, which sh ows a flux level of 
Fi.4GHz = (220+ 1 1) Jy dKuehretalJl98ll) . The hadronic radio 
surface profiles which are characterized by a smooth morphol- 
ogy fall short by roughly one order of magnitude even at impact 
paramete rs of some arc m inutes compared to observed profiles 
of lRottmann et alJ dl996l> . 

7. Predictions for next generation Cerenkov 
telescopes and GLAST 

There are three different scenarios predicting y-ray emission 
from objects like M 87, namely processed radiation of the rel- 
ativistic outflow, dark matter annihilation, and the hadronic 
scenario. The predictions of these types of model differ pre- 
dominantly in morphology, existing time-variability, and spec- 
tral signatures. The radio galaxy M 87 which is well within 
the field of view of the next generation Cerenkov telescopes 
MAGIC 1 , HESS 2 , and VERITAS 3 should therefore serve as 
a unique source for testing these scenarios. While the angu- 
lar resolution of these telescopes is comparable to the previ- 

1 http : //hegral . mppmu . mpg . de/MAGICWeb/ 

2 http : //www . mpi -hd . mpg . de/hfm/HESS/HESS . html 

3 http : //veritas . sao . arizona . edu/ 
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ously attained resolution, the flux sensitivities have strongly 
improved. These developments should allow for y-ray spec- 
troscopy by Cerenkov experiments in the near future, providing 
the opportunity of scrutinizing existing time-variation, and thus 
being able to constrain different y-ray emission scenarios. 

The LAT instrument onboard GLAST 4 will complement 
this research to even lower energies ranging from 20 MeV up 
to 300 GeV. Given a CRp population described by a single 
power-law spectral index a = ttSly extending from the GeV to 
TeV energy regime as well as a CRp scaling parameter Xqr p of 
Tabled we calculated the expected integrated y-ray flux above 
20 MeV. In the isobaric model, the y-ray flux estimates are 
T y {> 20 MeV)/(y ctrT 2 s" 1 ) = 6.0xl0" 8 , 1.3xl0" 7 , and2.3x 
10~ 7 for ajfy = 2.1, 2.2, and 2.27, respectively. This is well 
above the sensitivity limit of GLAST. The energy resolution 
of GLAST will even provide the possibility to disentangle the 
pion decay induced signature from inverse Compton emission 
of high-energetic electrons or positrons due to the energy reso- 
lution which is better than 10% and is sufficient to resolve the 
pion decay induced peak in the y-ray spectrum. 



8. Conclusion 

Using TeV y-ray detections of M 87 by the HEGRA collabo- 
ration, it is for the first time possible to constrain the CRp pop- 
ulation of an elliptical galaxy. By comparing to EGRET upper 
limits on the y-ray emission, we constrain the CRp spectral in- 
dex a<5.y, provided the y-ray emission is of hadronic origin and 
the population is described by a single power-law ranging from 
the GeV to TeV energy regime. 

By investigating three different models for the spatial dis- 
tribution of the CRp and applying those to realistic electron 
density and temperature profiles obtained from X-ray observa- 
tions, we calculate y-ray flux profiles resulting from hadronic 
CRp interaction s with the thermal am bient gas using an ana- 
lytic formalism jPfrommer & EnBlinl l2003). After convolving 
with the HEGRA point spread function, we compare the in- 
tegrated y-ray flux F y (E y > 730 GeV) for different annuli of 
equal solid angle elements centered on the source. Based on 
the available data we find good morphological agreement of all 
our models with these HEGRA Cerenkov observations. 

In the isobaric and adiabatic CRp model, the consequences 
for the CRp scaling parameter Xqr p drawn from normaliza- 
tion of our y-ray flux profiles to HEGRA observations yield 
slightly smaller values when comparing to our Galaxy with 
^CRp ~ 1) depending on the CRp spectral index. This is be- 
cause of the sensitivity of the observations to both the ISM of 
M 87 and the ICM of the central cooling flow region of Virgo, 
where Xcr p < 0.2 JPfrommer & EnBlin 2003). Especially for 
a GeV - 2.1 or lower momentum cutoffs of the CRp population 
due to Coulomb cooling processes, we obtain smaller contri- 
butions of CRp pressure to the ambient medium. By exploring 
our diffusion model and comparing our constraints on the CRp 
luminosity Lcr p to mechanical jet power estimates of M 87, 
we show the ability of TeV y-ray observations to constrain a 



combination of energy fraction of CRp escaping from the radio 
plasma and average jet power of the AGN. 

The expected radio emission by hadronically produced CRe 
is roughly one order of magnitude smaller compared to the syn- 
chrotron emission of the jet. Therefore it will be a challenge for 
future radio observations to disentangle the hadronic and jet 
emission components. Future Cerenkov observations should at 
least be able to severely constrain the parameter space of differ- 
ent y-ray emission scenarios. Finally, by investigating the y-ray 
flux in the energy regime of GLAST, we predict values which 
should allow to scrutinize this hadronic model in contrast to 
other scenarios providing us with the possibility of entering a 
new era of precision high energy cluster physics. 
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